!-----------------------------------------------------------------------
!
!  Date:  April 26, 2016 (Version 6)
!
!  Simple Physics Package
!
!  SIMPLE_PHYSICS includes large-scale precipitation, surface fluxes and
!  boundary-leyer mixing. The processes are time-split in that order.
!  A partially implicit formulation is used to foster numerical
!  stability. The routine assumes that the model levels are ordered
!  in a top-down approach, e.g. level 1 denotes the uppermost full model
!  level.
!
!  This routine is based on an implementation which was developed for
!  the NCAR Community Atmosphere Model (CAM). Adjustments for other
!  models may be necessary.
!
!  The routine provides both updates of the state variables u, v, T, q
!  (these are local copies of u,v,T,q within this physics routine) and
!  also collects their time tendencies. The latter might be used to
!  couple the physics and dynamics in a process-split way. For a
!  time-split coupling, the final state should be given to the
!  dynamical core for the next time step.
!
! Test:      0 = Reed and Jablonowski (2011) tropical cyclone test
!            1 = Moist baroclinic instability test
! RJ2012_precip:
!         true  = Turn on Reed and Jablonowski (2012) precip scheme
!         false = Turn off Reed and Jablonowski (2012) precip scheme
! TC_PBL_mod:
!         true  = Turn on George Bryan PBL mod for tropical cyclone test
!         false = Turn off George Bryan PBL mod (i.e., run as in Reed and Jablonowski (2012))
!
!  SUBROUTINE SIMPLE_PHYSICS(pcols, pver, dtime, lat, t, q, u, v, pmid, pint, pdel, ps, precl, test, RJ2012_precip, TC_PBL_mod)
!
!  Input variables:
!     pcols  - number of atmospheric columns (#)
!     pver   - number of model levels (#)
!     dtime  - time step (s)
!     lat    - latitude (radians)
!     t      - temperature at model levels (K)
!     q      - specific humidity at model levels (gm/gm)
!     u      - zonal wind at model levels (m/s)
!     v      - meridional wind at model levels (m/s)
!     pmid   - pressure at model levels (Pa)
!     pint   - pressure at interfaces (Pa)
!     pdel   - layer thickness (Pa)
!     ps     - surface pressure (Pa)
!     test   - test case to use for sea-surface temperatures
!     RJ2012_precip - RJ2012 precip flag
!     TC_PBL_mod    - PCL modification for TC test
!
!  Output variables:
!     Increments are added into t, q, u, v, pmid, pint, pdel, and ps
!     which are returned to the routine from which SIMPLE_PHYSICS was
!     called.  Precpitation is returned via precl.
!
!  Change log:
!  v2: removal of some NCAR CAM-specific 'use' associations
!  v3: corrected precl(i) computation, the precipitation rate is now
!      computed via a vertical integral, the previous single-level
!      computation in v2 was a bug
!  v3: corrected dtdt(i,1) computation, the term '-(i,1)' was missing
!      the temperature variable: '-t(i,1)'
!  v4: modified and enhanced parameter list to make the routine truly
!      standalone, the number of columns and vertical levels have been
!      added: pcols, pver
!  v4: 'ncol' has been removed, 'pcols' is used instead
!  v5: the sea surface temperature (SST) field Tsurf is now an array,
!      the SST now depends on the latitude
!  v5: addition of the latitude array 'lat' and the flag 'test' in the
!      parameter list
!      if test = 0: constant SST is used, correct setting for the
!                   tropical cyclone test
!      if test = 1: newly added latitude-dependent SST is used,
!                   correct setting for the moist baroclinic wave test
!                   with simple-physics
!  v6: addition of flags for a modified PBL for the TC test and
!      to turn off large-scale condensation scheme when using Kessler physics
!      Included virtual temperature in density calculation in PBL scheme
!      Also, included the virtual temperature, instead of temperature, for
!      the calculation of rho in the PBL scheme
!      (v6_1) Minor specification and generalization fixes.
!
! Reference: Reed, K. A. and C. Jablonowski (2012), Idealized tropical cyclone
!            simulations of intermediate complexity: A test case for AGCMs,
!            J. Adv. Model. Earth Syst., Vol. 4, M04001, doi:10.1029/2011MS000099
!-----------------------------------------------------------------------

subroutine simple_physics(pcols, pver, dtime, lat, t, q, u, v, pmid, pint, pdel, ps, &
                          dudt, dvdt, dtdt, dqdt, &
                          precl, test, RJ2012_precip, TC_PBL_mod)

  implicit none

  integer, parameter :: r8 = selected_real_kind(12)

  integer , intent(in) :: pcols                 ! Set number of atmospheric columns
  integer , intent(in) :: pver                  ! Set number of model levels
  real(r8), intent(in) :: dtime                 ! Set model physics timestep
  real(r8), intent(in) :: lat(pcols)            ! Latitude
  integer , intent(in) :: test                  ! Test number
  logical , intent(in) :: RJ2012_precip
  logical , intent(in) :: TC_PBL_mod

  real(r8), intent(inout) :: t   (pcols,pver)   ! Temperature at full-model level (K)
  real(r8), intent(inout) :: q   (pcols,pver)   ! Specific humidity at full-model level (kg/kg)
  real(r8), intent(inout) :: u   (pcols,pver)   ! Zonal wind at full-model level (m/s)
  real(r8), intent(inout) :: v   (pcols,pver)   ! Meridional wind at full-model level (m/s)
  real(r8), intent(inout) :: pmid(pcols,pver)   ! Pressure is full-model level (Pa)
  real(r8), intent(inout) :: pint(pcols,pver+1) ! Pressure at model interfaces (Pa)
  real(r8), intent(inout) :: pdel(pcols,pver)   ! Layer thickness (Pa)
  real(r8), intent(inout) :: ps  (pcols)        ! Surface pressue (Pa)

  real(r8), intent(out) :: precl(pcols)         ! Precipitation rate (m_water / s)

  integer i, k

  real(r8) gravit                               ! Gravity
  real(r8) rair                                 ! Gas constant for dry air
  real(r8) cpair                                ! Specific heat of dry air
  real(r8) latvap                               ! Latent heat of vaporization
  real(r8) rh2o                                 ! Gas constant for water vapor
  real(r8) epsilo                               ! Ratio of gas constant for dry air to that for vapor
  real(r8) zvir                                 ! Constant for virtual temp. calc. =(rh2o/rair) - 1
  real(r8) a                                    ! Reference Earth's Radius (m)
  real(r8) omega                                ! Reference rotation rate of the Earth (s^-1)
  real(r8) pi                                   ! pi

  ! Simple Physics Specific Constants
  real(r8) Tsurf(pcols)                         ! Sea surface temperature (constant for tropical cyclone)
                                                ! Tsurf needs to be dependent on latitude for the
                                                ! moist baroclinic wave test, adjust

  real(r8) SST_TC                               ! Sea surface temperature for tropical cyclone test
  real(r8) T0                                   ! Control temp for calculation of qsat
  real(r8) e0                                   ! Saturation vapor pressure at T0 for calculation of qsat
  real(r8) rhow                                 ! Density of Liquid Water
  real(r8) p0                                   ! Constant for calculation of potential temperature
  real(r8) Cd0                                  ! Constant for calculating Cd from Smith and Vogl 2008
  real(r8) Cd1                                  ! Constant for calculating Cd from Smith and Vogl 2008
  real(r8) Cm                                   ! Constant for calculating Cd from Smith and Vogl 2008
  real(r8) v20                                  ! Threshold wind speed for calculating Cd from Smith and Vogl 2008
  real(r8) C                                    ! Drag coefficient for sensible heat and evaporation
  real(r8) T00                                  ! Horizontal mean T at surface for moist baro test
  real(r8) u0                                   ! Zonal wind constant for moist baro test
  real(r8) latw                                 ! Halfwidth for baro test
  real(r8) eta0                                 ! Center of jets (hybrid) for baro test
  real(r8) etav                                 ! Auxiliary variable for baro test
  real(r8) q0                                   ! Maximum specific humidity for baro test
  real(r8) kappa                                ! von Karman constant

  real(r8), intent(inout) :: dtdt(pcols,pver)   ! Temperature tendency
  real(r8), intent(inout) :: dqdt(pcols,pver)   ! Specific humidity tendency
  real(r8), intent(inout) :: dudt(pcols,pver)   ! Zonal wind tendency
  real(r8), intent(inout) :: dvdt(pcols,pver)   ! Meridional wind tendency

  real(r8) tmp                                  ! Temporary
  real(r8) qsat                                 ! Saturation vapor pressure
  real(r8) qsats                                ! Saturation vapor pressure of SST

  ! Variables for Boundary Layer Calculation
  real(r8) wind(pcols       )                   ! Magnitude of Wind
  real(r8) Cd  (pcols       )                   ! Drag coefficient for momentum
  real(r8) Km  (pcols,pver+1)                   ! Eddy diffusivity for boundary layer calculations
  real(r8) Ke  (pcols,pver+1)                   ! Eddy diffusivity for boundary layer calculations
  real(r8) rho                                  ! Density at lower/upper interface
  real(r8) za  (pcols       )                   ! Heights at midpoints of first model level
  real(r8) zi  (pcols,pver+1)                   ! Heights at model interfaces
  real(r8) dlnpint                              ! Used for calculation of heights
  real(r8) pbltop                               ! Top of boundary layer
  real(r8) zpbltop                              ! Top of boundary layer for George Bryan Modifcation
  real(r8) pblconst                             ! Constant for the calculation of the decay of diffusivity
  real(r8) CA  (pcols,pver  )                   ! Matrix coefficents for PBL Scheme
  real(r8) CC  (pcols,pver  )                   ! Matrix coefficents for PBL Scheme
  real(r8) CE  (pcols,pver+1)                   ! Matrix coefficents for PBL Scheme
  real(r8) CAm (pcols,pver  )                   ! Matrix coefficents for PBL Scheme
  real(r8) CCm (pcols,pver  )                   ! Matrix coefficents for PBL Scheme
  real(r8) CEm (pcols,pver+1)                   ! Matrix coefficents for PBL Scheme
  real(r8) CFu (pcols,pver+1)                   ! Matrix coefficents for PBL Scheme
  real(r8) CFv (pcols,pver+1)                   ! Matrix coefficents for PBL Scheme
  real(r8) CFt (pcols,pver+1)                   ! Matrix coefficents for PBL Scheme
  real(r8) CFq (pcols,pver+1)                   ! Matrix coefficents for PBL Scheme


  ! Variable for Dry Mass Adjustment, this dry air adjustment is necessary to
  ! conserve the mass of the dry air

  gravit   = 9.80616_r8                         ! Gravity (9.80616 m/s^2)
  rair     = 287.0_r8                           ! Gas constant for dry air: 287 J/(kg K)
  cpair    = 1.0045e3_r8                        ! Specific heat of dry air: here we use 1004.5 J/(kg K)
  latvap   = 2.5e6_r8                           ! Latent heat of vaporization (J/kg)
  rh2o     = 461.5_r8                           ! Gas constant for water vapor: 461.5 J/(kg K)
  epsilo   = rair/rh2o                          ! Ratio of gas constant for dry air to that for vapor
  zvir     = (rh2o / rair) - 1                  ! Constant for virtual temp. calc. =(rh2o/rair) - 1 is approx. 0.608
  a        = 6371220.0_r8                       ! Reference Earth's Radius (m)
  omega    = 7.29212d-5                         ! Reference rotation rate of the Earth (s^-1)
  pi       = 4 * atan(1.0_r8)

  C        = 0.0011_r8                          ! From Smith and Vogl 2008
  SST_TC   = 302.15_r8                          ! Constant Value for SST
  T0       = 273.16_r8                          ! control temp for calculation of qsat
  e0       = 610.78_r8                          ! saturation vapor pressure at T0 for calculation of qsat
  rhow     = 1000.0_r8                          ! Density of Liquid Water
  Cd0      = 0.0007_r8                          ! Constant for Cd calc. Smith and Vogl 2008
  Cd1      = 0.000065_r8                        ! Constant for Cd calc. Smith and Vogl 2008
  Cm       = 0.002_r8                           ! Constant for Cd calc. Smith and Vogl 2008
  v20      = 20.0_r8                            ! Threshold wind speed for calculating Cd from Smith and Vogl 2008
  p0       = 100000.0_r8                        ! Constant for potential temp calculation
  pbltop   = 85000._r8                          ! Top of boundary layer in p
  zpbltop  = 1000._r8                           ! Top of boundary layer in z
  pblconst = 10000._r8                          ! Constant for the calculation of the decay of diffusivity
  T00      = 288.0_r8                           ! Horizontal mean T at surface for moist baro test
  u0       = 35.0_r8                            ! Zonal wind constant for moist baro test
  latw     = 2 * pi / 9.0_r8                    ! Halfwidth for baro test
  eta0     = 0.252_r8                           ! Center of jets (hybrid) for baro test
  etav     = (1 - eta0) * 0.5_r8 * pi           ! Auxiliary variable for baro test
  q0       = 0.021_r8                           ! Maximum specific humidity for baro test
  kappa    = 0.4_r8                             ! von Karman constant

  ! Calculate hydrostatic height
  do i = 1, pcols
    ! NOTE: ps(i) is identical to pint(i,pver+1), this is the correct sign (corrects typo in JAMES paper)
    dlnpint      = log(ps(i)) - log(pint(i,pver))
    za(i)        = rair / gravit * t(i,pver) * (1 + zvir * q(i,pver)) * 0.5_r8 * dlnpint
    zi(i,pver+1) = 0
  end do

  ! Set Sea Surface Temperature (constant for tropical cyclone)
  ! Tsurf needs to be dependent on latitude for moist baroclinic wave test
  ! Tsurf needs to be constant for tropical cyclone test
  select case (test)
  case (1) ! Moist Baroclinic Wave Test
    do i = 1, pcols
      Tsurf(i) = (T00 + pi*u0/rair * 1.5_r8 * sin(etav) * (cos(etav))**0.5_r8 *                 &
                 ((-2._r8*(sin(lat(i)))**6 * ((cos(lat(i)))**2 + 1._r8/3._r8) + 10._r8/63._r8)* &
                 u0 * (cos(etav))**1.5_r8  +                                                    &
                 (8._r8/5._r8*(cos(lat(i)))**3 * ((sin(lat(i)))**2 + 2._r8/3._r8) - pi/4._r8)*a*omega*0.5_r8 ))/ &
                 (1._r8+zvir*q0*exp(-(lat(i)/latw)**4))
    end do
  case (0) ! Tropical Cyclone Test
    do i = 1, pcols
      Tsurf(i) = SST_TC
    end do
  end select

  ! Set initial physics time tendencies and precipitation field to zero
  dtdt (:pcols,:pver) = 0 ! Initialize temperature tendency with zero
  dqdt (:pcols,:pver) = 0 ! Initialize specific humidity tendency with zero
  dudt (:pcols,:pver) = 0 ! Initialize zonal wind tendency with zero
  dvdt (:pcols,:pver) = 0 ! Initialize meridional wind tendency with zero
  precl(:pcols      ) = 0 ! Initialize precipitation rate with zero

  ! Large-Scale Condensation and Precipitation
  if (RJ2012_precip) then
    do k = 1, pver
      do i = 1, pcols
        qsat = epsilo * e0 / pmid(i,k) * exp(-latvap / rh2o * ((1.0_r8 / t(i,k)) - 1.0_r8 / T0))
        if (q(i,k) > qsat) then
          tmp  = 1.0_r8 / dtime * (q(i,k) - qsat) / (1 + (latvap / cpair) * (epsilo * latvap * qsat / (rair * t(i,k)**2)))
          dtdt (i,k) = dtdt(i,k) + latvap / cpair * tmp
          dqdt (i,k) = dqdt(i,k) - tmp
          precl(i  ) = precl(i) + tmp * pdel(i,k) / (gravit * rhow)
        end if
      end do
    end do

    ! Update moisture and temperature fields from Larger-Scale Precipitation Scheme
    do k = 1, pver
      do i = 1, pcols
        t(i,k) = t(i,k) + dtdt(i,k) * dtime
        q(i,k) = q(i,k) + dqdt(i,k) * dtime
      end do
    end do
  end if

  ! Turbulent mixing coefficients for the PBL mixing of horizontal momentum,
  ! sensible heat and latent heat
  !
  ! We are using Simplified Ekman theory to compute the diffusion coefficients
  ! Kx for the boundary-layer mixing. The Kx values are calculated at each time step
  ! and in each column.

  ! Compute magnitude of the wind and drag coeffcients for turbulence scheme
  do i = 1, pcols
    wind(i) = sqrt(u(i,pver)**2 + v(i,pver)**2)
  end do
  do i = 1, pcols
    if (wind(i) < v20) then
      Cd(i) = Cd0 + Cd1 * wind(i)
    else
      Cd(i) = Cm
    end if
  end do

  if (TC_PBL_mod) then ! Bryan TC PBL Modification
    do k = pver, 1, -1
      do i = 1, pcols
        dlnpint = log(pint(i,k+1)) - log(pint(i,k))
        zi(i,k) = zi(i,k+1) + rair / gravit * t(i,k) * (1 + zvir * q(i,k)) * dlnpint
        if (zi(i,k) <= zpbltop) then
          Km(i,k) = kappa * sqrt(Cd(i)) * wind(i) * zi(i,k) * (1 - zi(i,k) / zpbltop) * (1 - zi(i,k) / zpbltop)
          Ke(i,k) = kappa * sqrt(C)     * wind(i) * zi(i,k) * (1 - zi(i,k) / zpbltop) * (1 - zi(i,k) / zpbltop)
        else
          Km(i,k) = 0.0_r8
          Ke(i,k) = 0.0_r8
        end if
      end do
    end do
  else ! Reed and Jablonowski (2012) Configuration
    do k = 1, pver
      do i = 1, pcols
        if (pint(i,k) >= pbltop) then
          Km(i,k) = Cd(i) * wind(i) * za(i)
          Ke(i,k) = C * wind(i) * za(i)
        else
          Km(i,k) = Cd(i) * wind(i) * za(i) * exp(-(pbltop - pint(i,k))**2 / pblconst**2)
          Ke(i,k) = C     * wind(i) * za(i) * exp(-(pbltop - pint(i,k))**2 / pblconst**2)
        end if
      end do
    end do
  end if

  ! Update the state variables u, v, t, q with the surface fluxes at the
  ! lowest model level, this is done with an implicit approach
  ! see Reed and Jablonowski (JAMES, 2012)
  !
  ! Sea Surface Temperature Tsurf is constant for tropical cyclone test 5-1
  ! Tsurf needs to be dependent on latitude for the
  ! moist baroclinic wave test
  do i = 1, pcols
    qsats        = epsilo*e0/ps(i)*exp(-latvap/rh2o*((1._r8/Tsurf(i))-1._r8/T0))
    dudt(i,pver) = dudt(i,pver) + (u(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i))-u(i,pver))/dtime
    dvdt(i,pver) = dvdt(i,pver) + (v(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i))-v(i,pver))/dtime
    u   (i,pver) = u(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i))
    v   (i,pver) = v(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i))
    dtdt(i,pver) = dtdt(i,pver) +((t(i,pver)+C*wind(i)*Tsurf(i)*dtime/za(i))/(1._r8+C*wind(i)*dtime/za(i))-t(i,pver))/dtime
    t   (i,pver) = (t(i,pver)+C*wind(i)*Tsurf(i)*dtime/za(i))/(1._r8+C*wind(i)*dtime/za(i))
    tmp          = (q(i,pver) + C*wind(i) * qsats * dtime / za(i)) / (1 + C * wind(i) * dtime / za(i))
    dqdt(i,pver) = dqdt(i,pver) +(tmp - q(i,pver))/dtime
    q   (i,pver) = tmp ! Directly update q.
  end do

  ! Boundary layer mixing, see Reed and Jablonowski (JAMES, 2012)
  ! Calculate Diagonal Variables for Implicit PBL Scheme
  do k = 1, pver - 1
    do i = 1, pcols
      rho        = (pint(i,k+1)/(rair*(t(i,k+1)*(1._r8+zvir*q(i,k+1))+t(i,k)*(1._r8+zvir*q(i,k)))/2.0_r8))
      CAm(i,k)   = dtime/pdel(i,k)*gravit*gravit*Km(i,k+1)*rho*rho/(pmid(i,k+1)-pmid(i,k))
      CCm(i,k+1) = dtime/pdel(i,k+1)*gravit*gravit*Km(i,k+1)*rho*rho/(pmid(i,k+1)-pmid(i,k))
      CA (i,k)   = dtime/pdel(i,k)*gravit*gravit*Ke(i,k+1)*rho*rho/(pmid(i,k+1)-pmid(i,k))
      CC (i,k+1) = dtime/pdel(i,k+1)*gravit*gravit*Ke(i,k+1)*rho*rho/(pmid(i,k+1)-pmid(i,k))
    end do
  end do
  do i = 1, pcols
    CAm(i,pver)   = 0
    CCm(i,1)      = 0
    CEm(i,pver+1) = 0
    CA (i,pver)   = 0
    CC (i,1)      = 0
    CE (i,pver+1) = 0
    CFu(i,pver+1) = 0
    CFv(i,pver+1) = 0
    CFt(i,pver+1) = 0
    CFq(i,pver+1) = 0
  end do
  do i = 1, pcols
    do k = pver, 1, -1
      CE (i,k) = CC(i,k)/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1))
      CEm(i,k) = CCm(i,k)/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
      CFu(i,k) = (u(i,k)+CAm(i,k)*CFu(i,k+1))/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
      CFv(i,k) = (v(i,k)+CAm(i,k)*CFv(i,k+1))/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
      CFt(i,k) = ((p0/pmid(i,k))**(rair/cpair)*t(i,k)+CA(i,k)*CFt(i,k+1))/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1))
      CFq(i,k) = (q(i,k)+CA(i,k)*CFq(i,k+1))/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1))
    end do
  end do

  ! Calculate the updated temperaure and specific humidity and wind tendencies

  ! First we need to calculate the tendencies at the top model level
  do i = 1, pcols
    dudt(i,1) = dudt(i,1) + (CFu(i,1) - u(i,1)) / dtime
    dvdt(i,1) = dvdt(i,1) + (CFv(i,1) - v(i,1)) / dtime
    u   (i,1) = CFu (i,1)
    v   (i,1) = CFv (i,1)
    tmp       = CFt (i,1) * (pmid(i,1) / p0)**(rair / cpair)
    dtdt(i,1) = dtdt(i,1) + (tmp - t(i,1)) / dtime
    t   (i,1) = tmp
    dqdt(i,1) = dqdt(i,1) + (CFq (i,1) - q(i,1)) / dtime
    q   (i,1) = CFq(i,1)
  end do

  do i = 1, pcols
    do k = 2, pver
      dudt(i,k) = dudt(i,k) + (CEm(i,k) * u(i,k-1) + CFu(i,k) - u(i,k)) / dtime
      dvdt(i,k) = dvdt(i,k) + (CEm(i,k) * v(i,k-1) + CFv(i,k) - v(i,k)) / dtime
      u   (i,k) = CEm(i,k) * u(i,k-1) + CFu(i,k)
      v   (i,k) = CEm(i,k) * v(i,k-1) + CFv(i,k)
      tmp       = (CE(i,k) * t(i,k-1) * (p0 / pmid(i,k-1))**(rair / cpair) + CFt(i,k)) * (pmid(i,k) / p0)**(rair / cpair)
      dtdt(i,k) = dtdt(i,k) + (tmp - t(i,k)) / dtime
      t   (i,k) = tmp
      tmp       = CE(i,k) * q(i,k-1) + CFq(i,k)
      dqdt(i,k) = dqdt(i,k) + (tmp - q(i,k)) / dtime
      q   (i,k) = tmp
    end do
  end do

end subroutine simple_physics

